In [1]:
#Import necessary packages igraph, scvelo, scanpy, anndata, pandas, os, numpy
!pip install -U scvelo
from IPython.core.display import display, HTML
display(HTML("<style>.container { width:90% !important; }</style>"))
%matplotlib inline

import scvelo as scv
scv.logging.print_version()

#test whether igraph is imported
import igraph as ig
import matplotlib as plt
import loompy
import numpy as np
import scanpy as sc
import anndata
import pandas as pd
import os
import louvain
cwd = os.getcwd()
cwd
Requirement already up-to-date: scvelo in /usr/local/lib/python3.7/site-packages (0.1.25)
Requirement already satisfied, skipping upgrade: matplotlib>=2.2 in /usr/local/lib/python3.7/site-packages (from scvelo) (3.1.0)
Requirement already satisfied, skipping upgrade: scikit-learn>=0.21.2 in /usr/local/lib/python3.7/site-packages (from scvelo) (0.21.2)
Requirement already satisfied, skipping upgrade: numpy>=1.17 in /usr/local/lib/python3.7/site-packages (from scvelo) (1.18.3)
Requirement already satisfied, skipping upgrade: umap-learn>=0.3 in /usr/local/lib/python3.7/site-packages (from scvelo) (0.3.9)
Requirement already satisfied, skipping upgrade: pandas>=0.23 in /usr/local/lib/python3.7/site-packages (from scvelo) (0.24.2)
Requirement already satisfied, skipping upgrade: scanpy>=1.4 in /usr/local/lib/python3.7/site-packages (from scvelo) (1.4.3)
Requirement already satisfied, skipping upgrade: anndata>=0.6.21 in /usr/local/lib/python3.7/site-packages (from scvelo) (0.6.21)
Requirement already satisfied, skipping upgrade: scipy>=1.0 in /usr/local/lib/python3.7/site-packages (from scvelo) (1.4.1)
Requirement already satisfied, skipping upgrade: loompy>=2.0.12 in /usr/local/lib/python3.7/site-packages (from scvelo) (3.0.6)
Requirement already satisfied, skipping upgrade: cycler>=0.10 in /usr/local/lib/python3.7/site-packages (from matplotlib>=2.2->scvelo) (0.10.0)
Requirement already satisfied, skipping upgrade: pyparsing!=2.0.4,!=2.1.2,!=2.1.6,>=2.0.1 in /usr/local/lib/python3.7/site-packages (from matplotlib>=2.2->scvelo) (2.4.0)
Requirement already satisfied, skipping upgrade: python-dateutil>=2.1 in /usr/local/lib/python3.7/site-packages (from matplotlib>=2.2->scvelo) (2.8.0)
Requirement already satisfied, skipping upgrade: kiwisolver>=1.0.1 in /usr/local/lib/python3.7/site-packages (from matplotlib>=2.2->scvelo) (1.1.0)
Requirement already satisfied, skipping upgrade: joblib>=0.11 in /usr/local/lib/python3.7/site-packages (from scikit-learn>=0.21.2->scvelo) (0.13.2)
Requirement already satisfied, skipping upgrade: numba>=0.37 in /usr/local/lib/python3.7/site-packages (from umap-learn>=0.3->scvelo) (0.44.0)
Requirement already satisfied, skipping upgrade: pytz>=2011k in /usr/local/lib/python3.7/site-packages (from pandas>=0.23->scvelo) (2019.1)
Requirement already satisfied, skipping upgrade: natsort in /usr/local/lib/python3.7/site-packages (from scanpy>=1.4->scvelo) (6.0.0)
Requirement already satisfied, skipping upgrade: tables in /usr/local/lib/python3.7/site-packages (from scanpy>=1.4->scvelo) (3.5.2)
Requirement already satisfied, skipping upgrade: networkx in /usr/local/lib/python3.7/site-packages (from scanpy>=1.4->scvelo) (2.3)
Requirement already satisfied, skipping upgrade: patsy in /usr/local/lib/python3.7/site-packages (from scanpy>=1.4->scvelo) (0.5.1)
Requirement already satisfied, skipping upgrade: seaborn in /usr/local/lib/python3.7/site-packages (from scanpy>=1.4->scvelo) (0.9.0)
Requirement already satisfied, skipping upgrade: statsmodels in /usr/local/lib/python3.7/site-packages (from scanpy>=1.4->scvelo) (0.9.0)
Requirement already satisfied, skipping upgrade: h5py in /usr/local/lib/python3.7/site-packages (from scanpy>=1.4->scvelo) (2.9.0)
Requirement already satisfied, skipping upgrade: tqdm in /usr/local/lib/python3.7/site-packages (from scanpy>=1.4->scvelo) (4.32.1)
Requirement already satisfied, skipping upgrade: numpy-groupies in /usr/local/lib/python3.7/site-packages (from loompy>=2.0.12->scvelo) (0+unknown)
Requirement already satisfied, skipping upgrade: click in /usr/local/lib/python3.7/site-packages (from loompy>=2.0.12->scvelo) (7.0)
Requirement already satisfied, skipping upgrade: setuptools in /usr/local/lib/python3.7/site-packages (from loompy>=2.0.12->scvelo) (46.0.0)
Requirement already satisfied, skipping upgrade: six in /usr/local/lib/python3.7/site-packages (from cycler>=0.10->matplotlib>=2.2->scvelo) (1.12.0)
Requirement already satisfied, skipping upgrade: llvmlite>=0.29.0 in /usr/local/lib/python3.7/site-packages (from numba>=0.37->umap-learn>=0.3->scvelo) (0.29.0)
Requirement already satisfied, skipping upgrade: numexpr>=2.6.2 in /usr/local/lib/python3.7/site-packages (from tables->scanpy>=1.4->scvelo) (2.6.9)
Requirement already satisfied, skipping upgrade: mock>=2.0 in /usr/local/lib/python3.7/site-packages (from tables->scanpy>=1.4->scvelo) (3.0.5)
Requirement already satisfied, skipping upgrade: decorator>=4.3.0 in /usr/local/lib/python3.7/site-packages (from networkx->scanpy>=1.4->scvelo) (4.4.0)
Running scvelo 0.1.25 (python 3.7.7) on 2020-05-09 17:41.
Out[1]:
'/Users/whippoorwill/Desktop/Jupyter'
In [72]:
# from scvelo.tools.velocity_embedding import velocity_embedding
# from scvelo.tools.utils import groups_to_bool
# from scvelo.plotting.utils import default_basis, default_size, default_color, get_components, savefig_or_show, make_unique_list,get_basis, velocity_embedding_changed, get_ax
# from scvelo.pl.velocity_embedding_grid import compute_velocity_on_grid
# from scvelo.pl.scatter import scatter
# from scvelo.pl.docs import doc_scatter, doc_params

# from matplotlib import rcParams
# import matplotlib.pyplot as pl
# import numpy as np


# @doc_params(scatter=doc_scatter)
# def velocity_embedding_stream(adata, basis=None, vkey='velocity', density=None, smooth=None, min_mass=None, cutoff_perc=None,
#                               arrow_color=None, linewidth=None, n_neighbors=None, recompute=None, color=None, use_raw=None,
#                               layer=None, color_map=None, colorbar=True, palette=None, size=None, alpha=.3, perc=None,
#                               X=None, V=None, X_grid=None, V_grid=None, sort_order=True, groups=None, components=None,
#                               legend_loc='on data', legend_fontsize=None, legend_fontweight=None, xlabel=None,
#                               ylabel=None, title=None, fontsize=None, figsize=None, dpi=None, frameon=None, show=True,
#                               save=None, ax=None, ncols=None, **kwargs):
#     """\
#     Stream plot of velocities on the embedding.
#     Arguments
#     ---------
#     adata: :class:`~anndata.AnnData`
#         Annotated data matrix.
#     density: `float` (default: 1)
#         Amount of velocities to show - 0 none to 1 all
#     smooth: `float` (default: 0.5)
#         Multiplication factor for scale in Gaussian kernel around grid point.
#     min_mass: `float` (default: 1)
#         Minimum threshold for mass to be shown. It can range between 0 (all velocities) and 5 (large velocities only).
#     cutoff_perc: `float` (default: `None`)
#         If set, mask small velocities below a percentile threshold (range between 0 and 100).
#     linewidth: `float` (default: 1)
#         Line width for streamplot.
#     n_neighbors: `int` (default: None)
#         Number of neighbors to consider around grid point.
#     X: `np.ndarray` (default: None)
#         Embedding grid point coordinates
#     V: `np.ndarray` (default: None)
#         Embedding grid velocity coordinates
#     {scatter}
#     Returns
#     -------
#         `matplotlib.Axis` if `show==False`
#     """
#     basis = default_basis(adata, **kwargs) if basis is None else get_basis(adata, basis)
#     vkey = [key for key in list(adata.layers.keys()) if 'velocity' in key and '_u' not in key] if vkey == 'all' else vkey
#     color, color_map = kwargs.pop('c', color), kwargs.pop('cmap', color_map)
#     colors, layers, vkeys = make_unique_list(color, allow_array=True), make_unique_list(layer), make_unique_list(vkey)

#     if V is None:
#         for key in vkeys:
#             if recompute or velocity_embedding_changed(adata, basis=basis, vkey=key):
#                 velocity_embedding(adata, basis=basis, vkey=key)

#     color, layer, vkey = colors[0], layers[0], vkeys[0]
#     color = default_color(adata) if color is None else color

#     if X_grid is None or V_grid is None:
#         _adata = adata[groups_to_bool(adata, groups, groupby=color)] \
#             if groups is not None and color in adata.obs.keys() else adata
#         X_emb = np.array(_adata.obsm['X_' + basis][:, get_components(components, basis)]) if X is None else X[:, :2]
#         V_emb = np.array(_adata.obsm[vkey + '_' + basis][:, get_components(components, basis)]) if V is None else V[:, :2]
#         X_grid, V_grid = compute_velocity_on_grid(X_emb=X_emb, V_emb=V_emb, density=1, smooth=smooth, min_mass=min_mass,
#                                                   n_neighbors=n_neighbors, autoscale=False, adjust_for_stream=True,
#                                                   cutoff_perc=cutoff_perc)
#         lengths = np.sqrt((V_grid ** 2).sum(0))
#         linewidth = 1 if linewidth is None else linewidth
#         linewidth *= 2 * lengths / lengths[~np.isnan(lengths)].max()

#     scatter_kwargs = {"basis": basis, "perc": perc, "use_raw": use_raw, "sort_order": sort_order, "alpha": alpha,
#                       "components": components, "legend_loc": legend_loc, "groups": groups,
#                       "legend_fontsize": legend_fontsize, "legend_fontweight": legend_fontweight, "palette": palette,
#                       "color_map": color_map, "frameon": frameon, "xlabel": xlabel, "ylabel": ylabel,
#                       "colorbar": colorbar, "dpi": dpi, "fontsize": fontsize, "show": False, "save": False}

#     multikey = colors if len(colors) > 1 else layers if len(layers) > 1 else vkeys if len(vkeys) > 1 else None
#     if multikey is not None:
#         if title is None: title = list(multikey)
#         elif isinstance(title, (list, tuple)): title *= int(np.ceil(len(multikey) / len(title)))
#         ncols = len(multikey) if ncols is None else min(len(multikey), ncols)
#         nrows = int(np.ceil(len(multikey) / ncols))
#         figsize = rcParams['figure.figsize'] if figsize is None else figsize
#         ax = []
#         for i, gs in enumerate(
#                 pl.GridSpec(nrows, ncols, pl.figure(None, (figsize[0] * ncols, figsize[1] * nrows), dpi=dpi))):
#             if i < len(multikey):
#                 ax.append(velocity_embedding_stream(adata, density=density, size=size, smooth=smooth, n_neighbors=n_neighbors,
#                                                     linewidth=linewidth, ax=pl.subplot(gs),
#                                                     color=colors[i] if len(colors) > 1 else color,
#                                                     layer=layers[i] if len(layers) > 1 else layer,
#                                                     vkey=vkeys[i] if len(vkeys) > 1 else vkey,
#                                                     title=title[i] if isinstance(title, (list, tuple)) else title,
#                                                     X_grid=None if len(vkeys) > 1 else X_grid,
#                                                     V_grid=None if len(vkeys) > 1 else V_grid, **scatter_kwargs, **kwargs))
#         savefig_or_show(dpi=dpi, save=save, show=show)
#         if not show: return ax

#     else:
#         ax, show = get_ax(ax, show, figsize, dpi)
#         density = 1 if density is None else density
#         stream_kwargs = {"linewidth": linewidth, "density": 2 * density, "zorder": 3,
#                          "color": 'k' if arrow_color is None else arrow_color}
#         for arg in list(kwargs):
#             if arg in stream_kwargs: stream_kwargs.update({arg: kwargs[arg]})
#             else: scatter_kwargs.update({arg: kwargs[arg]})

#         ax.streamplot(X_grid[0], X_grid[1], V_grid[0], V_grid[1], **stream_kwargs)

#         size = 8 * default_size(adata) if size is None else size
#         ax = scatter(adata, layer=layer, color=color, size=size, title=title, ax=ax, zorder=0, **scatter_kwargs)

#         savefig_or_show(dpi=dpi, save=save, show=show)
#         if not show: return ax
---------------------------------------------------------------------------
ImportError                               Traceback (most recent call last)
<ipython-input-72-baa4132f185e> in <module>
      1 from scvelo.tools.velocity_embedding import velocity_embedding
      2 from scvelo.tools.utils import groups_to_bool
----> 3 from scvelo.plotting.utils import default_basis, default_size, default_color, get_components, savefig_or_show, make_unique_list,get_basis, velocity_embedding_changed, get_ax
      4 from scvelo.pl.velocity_embedding_grid import compute_velocity_on_grid
      5 from scvelo.pl.scatter import scatter

ImportError: cannot import name 'get_basis' from 'scvelo.plotting.utils' (/usr/local/lib/python3.7/site-packages/scvelo/plotting/utils.py)

Import the data file from inside the Jupyter folder as "adata"

In [2]:
adata = scv.read("mgAVM02/adata_mgAVM02.h5ad",cache = True)

If you get a message that "variable names are not unique", run the first line (make_unique). Also set verbosity and figure params for good visualization.

In [3]:
#adata.var_names_make_unique()
scv.settings.verbosity = 4  # show errors(0), warnings(1), info(2) and hints(3)
scv.settings.set_figure_params('scvelo')  # for beautified visualization
In [4]:
# show proportions of spliced/unspliced abundances
scv.utils.show_proportions(adata)
adata
Abundance of ['spliced', 'unspliced', 'ambiguous']: [0.5  0.41 0.09]
Out[4]:
AnnData object with n_obs × n_vars = 12330 × 5976 
    obs: 'nCount_RNA', 'percent.mito', 'age', 'condition', 'finalclusters', 'sample_description', 'sex', 'initial_size_spliced', 'initial_size_unspliced', 'initial_size', 'n_counts'
    var: 'Accession', 'Chromosome', 'End', 'Start', 'Strand', 'means', 'dispersions', 'dispersions_norm', 'velocity_offset', 'velocity_offset2', 'velocity_gamma', 'velocity_r2', 'velocity_genes', 'velocity_score'
    uns: 'age_colors', 'condition_colors', 'finalclusters_colors', 'neighbors', 'rank_velocity_genes', 'velocity_graph', 'velocity_graph_neg'
    obsm: 'X_umap', 'X_pca', 'velocity_umap'
    layers: 'Ms', 'Mu', 'ambiguous', 'matrix', 'spliced', 'unspliced', 'variance_velocity', 'velocity'
In [5]:
print (adata.obs_keys) #find out whether there is any annotation yet in adata
<bound method AnnData.obs_keys of AnnData object with n_obs × n_vars = 12330 × 5976 
    obs: 'nCount_RNA', 'percent.mito', 'age', 'condition', 'finalclusters', 'sample_description', 'sex', 'initial_size_spliced', 'initial_size_unspliced', 'initial_size', 'n_counts'
    var: 'Accession', 'Chromosome', 'End', 'Start', 'Strand', 'means', 'dispersions', 'dispersions_norm', 'velocity_offset', 'velocity_offset2', 'velocity_gamma', 'velocity_r2', 'velocity_genes', 'velocity_score'
    uns: 'age_colors', 'condition_colors', 'finalclusters_colors', 'neighbors', 'rank_velocity_genes', 'velocity_graph', 'velocity_graph_neg'
    obsm: 'X_umap', 'X_pca', 'velocity_umap'
    layers: 'Ms', 'Mu', 'ambiguous', 'matrix', 'spliced', 'unspliced', 'variance_velocity', 'velocity'>
In [6]:
print(max(adata.var['dispersions_norm']))
print(min(adata.var['dispersions_norm']))
print(sum(adata.var['dispersions_norm'])/len(adata.var['dispersions_norm']))
len(adata.var['dispersions_norm'])
x = sorted(adata.var['dispersions_norm'],reverse=True)
22.81192398071289
-0.0674969032406807
0.6650465992924447
In [7]:
adata
adata.obs['condition'][1:3]
Out[7]:
index
LD5RC:AAACCCAGTGCCTTTCx    Control
LD5RC:AAACGAAAGGTGGGTTx    Control
Name: condition, dtype: category
Categories (2, object): [Control, Deprived]
In [8]:
a = adata.var['velocity_genes']
b = adata.var['velocity_gamma']
c = adata.var['velocity_r2']
d = adata.layers['velocity']
In [9]:
adata.layers['velocity']
sum(adata.var['velocity_genes'])
Out[9]:
560

If you want to look at the velocity of a single gene, put it in the code below.

In [10]:
scv.pl.velocity(adata,var_names = "Cx3cr1",show=True,color='condition')
scv.pl.velocity(adata,var_names = "P2ry12",show=True,color = 'age')
scv.pl.velocity(adata,var_names = "Nrxn1",show=True,color = 'finalclusters')
scv.pl.velocity(adata,var_names = "Gria2",show=True,color = 'finalclusters')
scv.pl.velocity(adata[adata.obs["age"] == "P5",:],var_names = "P2ry12",show = True, color = 'finalclusters')
scv.pl.velocity(adata[adata.obs["age"] == "P7",:],var_names = "P2ry12",show = True, color = 'finalclusters')
Trying to set attribute `.obs` of view, making a copy.
Trying to set attribute `.obs` of view, making a copy.
Trying to set attribute `.obs` of view, making a copy.
Trying to set attribute `.obs` of view, making a copy.
Trying to set attribute `.obs` of view, making a copy.
Trying to set attribute `.obs` of view, making a copy.
Trying to set attribute `.obs` of view, making a copy.
Trying to set attribute `.obs` of view, making a copy.
In [11]:
bdata = adata[adata.obs['finalclusters'] == '3',:]

scv.pl.velocity(bdata,var_names = "Mpeg1",show = True, color = 'condition')
scv.pl.velocity(bdata,var_names = "Nrxn1",show = True, color = 'condition')
scv.pl.velocity(bdata,var_names = "Shank1",show = True, color = 'condition')
scv.pl.velocity(bdata,var_names = "Gria2",show = True, color = 'condition')
scv.pl.velocity(bdata,var_names = "Meg3",show = True, color = 'condition')
scv.pl.velocity(bdata,var_names = "Map1b",show = True, color = 'condition')
scv.pl.velocity(bdata,var_names = "Mast1",show = True, color = 'condition')
scv.pl.velocity(bdata,var_names = "Sparcl1",show = True, color = 'condition')
scv.pl.velocity(bdata,var_names = "Lamp1",show = True, color = 'condition')
Trying to set attribute `.obs` of view, making a copy.
Trying to set attribute `.obs` of view, making a copy.
Trying to set attribute `.obs` of view, making a copy.
Trying to set attribute `.obs` of view, making a copy.
In [12]:
sc.pl.umap(adata,color='condition',palette=['green','pink']) #show umap plot
In [13]:
sc.pl.umap(adata,color = 'age',palette=['blue','green']) #show umap
In [89]:
#adata.uns
In [90]:
#adata.var['dispersions_norm']
In [120]:
from functools import partial
def _wraps_plot(wrapper, func):
    annots_orig = {k: v for k, v in wrapper.__annotations__.items() if k not in {'self', 'kwargs'}}  # 'adata',
    annots = {k: v for k, v in func.__annotations__.items()}
    wrapper.__annotations__ = {**annots, **annots_orig}
    wrapper.__wrapped__ = func
    return wrapper


_wraps_plot_scatter = partial(_wraps_plot, func=scatter)
_wraps_plot_velocity_embedding = partial(_wraps_plot, func=velocity_embedding)
_wraps_plot_velocity_embedding_grid = partial(_wraps_plot, func=velocity_embedding_grid)
_wraps_plot_velocity_embedding_stream = partial(_wraps_plot, func=velocity_embedding_stream)
_wraps_plot_velocity_graph = partial(_wraps_plot, func=velocity_graph)
_wraps_plot_hist = partial(_wraps_plot, func=hist)

def gridspec(ncols=4, nrows=1, figsize=None, dpi=None):
    if figsize is None: figsize = rcParams['figure.figsize']
    gs = GridSpec(nrows, ncols, pl.figure(None, (figsize[0] * ncols, figsize[1] * nrows), dpi=dpi))
    return gs


class GridSpec:
    def __init__(self, ncols=4, nrows=1, figsize=None, dpi=None, **scatter_kwargs):
        """Specifies the geometry of the grid that a subplots can be placed in
        Example
        .. code:: python
            with scv.GridSpec() as pl:
                pl.scatter(adata, basis='pca')
                pl.scatter(adata, basis='umap')
                pl.hist(adata.obs.initial_size)
        Parameters
        ----------
        ncols: `int` (default: 4)
            Number of panels per row.
        nrows: `int` (default: 1)
            Number of panels per column.
        figsize: tuple (default: `None`)
            Figure size.
        dpi: `int` (default: `None`)
            Figure dpi.
        scatter_kwargs:
            Arguments to be passed to all scatter plots, e.g. `frameon=False`.
        """
        self.ncols, self.nrows, self.figsize, self.dpi = ncols, nrows, figsize, dpi
        self.scatter_kwargs = scatter_kwargs
        self.scatter_kwargs.update({'show': False})
        self.get_new_grid()
        self.new_row = None

    def __enter__(self):
        return self

    def __exit__(self, exc_type, exc_val, exc_tb):
        if self.new_row and self.count < self.max_count:
            ax = pl.subplot(self.gs[self.max_count - 1])
            ax.axis('off')
        pl.show()

    def get_new_grid(self):
        self.gs = gridspec(self.ncols, self.nrows, self.figsize, self.dpi)
        geo = self.gs[0].get_geometry()
        self.max_count, self.count, self.new_row = geo[0] * geo[1], 0, True

    def get_ax(self):
        if self.count >= self.max_count:
            self.get_new_grid()
        self.count += 1
        return pl.subplot(self.gs[self.count - 1])

    def get_kwargs(self, kwargs=None):
        _kwargs = self.scatter_kwargs.copy()
        if kwargs is not None:
            _kwargs.update(kwargs)
        _kwargs.update({'ax': self.get_ax(), 'show': False})
        return _kwargs

    @_wraps_plot_scatter
    def scatter(self, adata, **kwargs):
        return scatter(adata, **self.get_kwargs(kwargs))

    @_wraps_plot_velocity_embedding
    def velocity_embedding(self, adata, **kwargs):
        return velocity_embedding(adata, **self.get_kwargs(kwargs))

    @_wraps_plot_velocity_embedding_grid
    def velocity_embedding_grid(self, adata, **kwargs):
        return velocity_embedding_grid(adata, **self.get_kwargs(kwargs))

    @_wraps_plot_velocity_embedding_stream
    def velocity_embedding_stream(self, adata, **kwargs):
        return velocity_embedding_stream(adata, **self.get_kwargs(kwargs))

    @_wraps_plot_velocity_graph
    def velocity_embedding(self, adata, **kwargs):
        return velocity_embedding(adata, **self.get_kwargs(kwargs))

    @_wraps_plot_velocity_graph
    def velocity_graph(self, adata, **kwargs):
        return velocity_graph(adata, **self.get_kwargs(kwargs))

    @_wraps_plot_hist
    def hist(self, array, **kwargs):
        return hist(array, **self.get_kwargs(kwargs))
    
    
from sklearn.neighbors import NearestNeighbors
from scipy.stats import norm as normal
def compute_velocity_on_grid(X_emb, V_emb, density=None, smooth=None, n_neighbors=None, min_mass=None, autoscale=True,
                             adjust_for_stream=False, cutoff_perc=None):
    # remove invalid cells
    idx_valid = np.isfinite(X_emb.sum(1) + V_emb.sum(1))
    X_emb = X_emb[idx_valid]
    V_emb = V_emb[idx_valid]

    # prepare grid
    n_obs, n_dim = X_emb.shape
    density = 1 if density is None else density
    smooth = .5 if smooth is None else smooth

    grs = []
    for dim_i in range(n_dim):
        m, M = np.min(X_emb[:, dim_i]), np.max(X_emb[:, dim_i])
        m = m - .01 * np.abs(M - m)
        M = M + .01 * np.abs(M - m)
        gr = np.linspace(m, M, int(50 * density))
        grs.append(gr)

    meshes_tuple = np.meshgrid(*grs)
    X_grid = np.vstack([i.flat for i in meshes_tuple]).T

    # estimate grid velocities
    if n_neighbors is None: n_neighbors = int(n_obs/50)
    nn = NearestNeighbors(n_neighbors=n_neighbors, n_jobs=-1)
    nn.fit(X_emb)
    dists, neighs = nn.kneighbors(X_grid)

    scale = np.mean([(g[1] - g[0]) for g in grs]) * smooth
    weight = normal.pdf(x=dists, scale=scale)
    p_mass = weight.sum(1)

    V_grid = (V_emb[neighs] * weight[:, :, None]).sum(1) / np.maximum(1, p_mass)[:, None]
    if min_mass is None: min_mass = 1

    if adjust_for_stream:
        X_grid = np.stack([np.unique(X_grid[:, 0]), np.unique(X_grid[:, 1])])
        ns = int(np.sqrt(len(V_grid[:, 0])))
        V_grid = V_grid.T.reshape(2, ns, ns)

        mass = np.sqrt((V_grid ** 2).sum(0))
        min_mass = 10 ** (min_mass - 6)  # default min_mass = 1e-5
        min_mass = np.clip(min_mass, None, np.max(mass) * .9)
        cutoff = mass.reshape(V_grid[0].shape) < min_mass

        if cutoff_perc is None: cutoff_perc = 5
        length = np.sum(np.mean(np.abs(V_emb[neighs]), axis=1), axis=1).T.reshape(ns, ns)
        cutoff |= length < np.percentile(length, cutoff_perc)

        V_grid[0][cutoff] = np.nan
    else:
        min_mass *= np.percentile(p_mass, 99) / 100
        X_grid, V_grid = X_grid[p_mass > min_mass], V_grid[p_mass > min_mass]

        if autoscale: V_grid /= 3 * quiver_autoscale(X_grid, V_grid)

    return X_grid, V_grid

 
def velocity_embedding_changed(adata, basis, vkey):
    if 'X_' + basis not in adata.obsm.keys(): changed = False
    else:
        changed = vkey + '_' + basis not in adata.obsm_keys()
        if vkey + '_settings' in adata.uns.keys():
            sett = adata.uns[vkey + '_settings']
            changed |= 'embeddings' not in sett or basis not in sett['embeddings']
    return changed

def check_basis(adata, basis):
    if basis in adata.obsm.keys() and 'X_' + basis not in adata.obsm.keys():
        adata.obsm['X_' + basis] = adata.obsm[basis]
        logg.info('Renamed', '\'' + basis + '\'', 'to convention', '\'X_' + basis + '\' (adata.obsm).')


def get_basis(adata, basis):
    if isinstance(basis, str) and basis.startswith('X_'):
        basis = basis[2:]
    check_basis(adata, basis)
    return basis

def velocity_embedding_stream(adata, basis=None, vkey='velocity', density=None, smooth=None, min_mass=None, cutoff_perc=None,
                              arrow_color=None, linewidth=None, n_neighbors=None, recompute=None, color=None, use_raw=None,
                              layer=None, color_map=None, colorbar=True, palette=None, size=None, alpha=.3, perc=None,
                              X=None, V=None, X_grid=None, V_grid=None, sort_order=True, groups=None, components=None,
                              legend_loc='on data', legend_fontsize=None, legend_fontweight=None, xlabel=None,
                              ylabel=None, title=None, fontsize=None, figsize=None, dpi=None, frameon=None, show=True,
                              save=None, ax=None, ncols=None, **kwargs):
    """\
    Stream plot of velocities on the embedding.
    Arguments
    ---------
    adata: :class:`~anndata.AnnData`
        Annotated data matrix.
    density: `float` (default: 1)
        Amount of velocities to show - 0 none to 1 all
    smooth: `float` (default: 0.5)
        Multiplication factor for scale in Gaussian kernel around grid point.
    min_mass: `float` (default: 1)
        Minimum threshold for mass to be shown. It can range between 0 (all velocities) and 5 (large velocities only).
    cutoff_perc: `float` (default: `None`)
        If set, mask small velocities below a percentile threshold (range between 0 and 100).
    linewidth: `float` (default: 1)
        Line width for streamplot.
    n_neighbors: `int` (default: None)
        Number of neighbors to consider around grid point.
    X: `np.ndarray` (default: None)
        Embedding grid point coordinates
    V: `np.ndarray` (default: None)
        Embedding grid velocity coordinates
    {scatter}
    Returns
    -------
        `matplotlib.Axis` if `show==False`
    """
    
    basis = default_basis(adata, **kwargs) if basis is None else get_basis(adata, basis)
    vkey = [key for key in list(adata.layers.keys()) if 'velocity' in key and '_u' not in key] if vkey == 'all' else vkey
    color, color_map = kwargs.pop('c', color), kwargs.pop('cmap', color_map)
    colors, layers, vkeys = make_unique_list(color, allow_array=True), make_unique_list(layer), make_unique_list(vkey)

    if V is None:
        for key in vkeys:
            if recompute or velocity_embedding_changed(adata, basis=basis, vkey=key):
                velocity_embedding(adata, basis=basis, vkey=key)

    color, layer, vkey = colors[0], layers[0], vkeys[0]
    color = default_color(adata) if color is None else color

    if X_grid is None or V_grid is None:
        _adata = adata[groups_to_bool(adata, groups, groupby=color)] \
            if groups is not None and color in adata.obs.keys() else adata
        X_emb = np.array(_adata.obsm['X_' + basis][:, get_components(components, basis)]) if X is None else X[:, :2]
        V_emb = np.array(_adata.obsm[vkey + '_' + basis][:, get_components(components, basis)]) if V is None else V[:, :2]
        X_grid, V_grid = compute_velocity_on_grid(X_emb=X_emb, V_emb=V_emb, density=1, smooth=smooth, min_mass=min_mass,
                                                  n_neighbors=n_neighbors, autoscale=False, adjust_for_stream=True,
                                                  cutoff_perc=cutoff_perc)
        lengths = np.sqrt((V_grid ** 2).sum(0))
        linewidth = 1 if linewidth is None else linewidth
        linewidth *= 2 * lengths / lengths[~np.isnan(lengths)].max()

    scatter_kwargs = {"basis": basis, "perc": perc, "use_raw": use_raw, "sort_order": sort_order, "alpha": alpha,
                      "components": components, "legend_loc": legend_loc, "groups": groups,
                      "legend_fontsize": legend_fontsize, "legend_fontweight": legend_fontweight, "palette": palette,
                      "color_map": color_map, "frameon": frameon, "xlabel": xlabel, "ylabel": ylabel,
                      "colorbar": colorbar, "dpi": dpi, "fontsize": fontsize, "show": False, "save": False}

    multikey = colors if len(colors) > 1 else layers if len(layers) > 1 else vkeys if len(vkeys) > 1 else None
    if multikey is not None:
        if title is None: title = list(multikey)
        elif isinstance(title, (list, tuple)): title *= int(np.ceil(len(multikey) / len(title)))
        ncols = len(multikey) if ncols is None else min(len(multikey), ncols)
        nrows = int(np.ceil(len(multikey) / ncols))
        figsize = rcParams['figure.figsize'] if figsize is None else figsize
        ax = []
        for i, gs in enumerate(
                GridSpec(nrows, ncols, pl.figure(None, (figsize[0] * ncols, figsize[1] * nrows), dpi=dpi))):
            if i < len(multikey):
                ax.append(velocity_embedding_stream(adata, density=density, size=size, smooth=smooth, n_neighbors=n_neighbors,
                                                    linewidth=linewidth, ax=pl.subplot(gs),
                                                    color=colors[i] if len(colors) > 1 else color,
                                                    layer=layers[i] if len(layers) > 1 else layer,
                                                    vkey=vkeys[i] if len(vkeys) > 1 else vkey,
                                                    title=title[i] if isinstance(title, (list, tuple)) else title,
                                                    X_grid=None if len(vkeys) > 1 else X_grid,
                                                    V_grid=None if len(vkeys) > 1 else V_grid, **scatter_kwargs, **kwargs))
        savefig_or_show(dpi=dpi, save=save, show=show)
        if not show: return ax

    else:
        ax, show = get_ax(ax, show, figsize, dpi)
        density = 1 if density is None else density
        stream_kwargs = {"linewidth": linewidth, "density": 2 * density, "zorder": 3,
                         "color": 'k' if arrow_color is None else arrow_color}
        for arg in list(kwargs):
            if arg in stream_kwargs: stream_kwargs.update({arg: kwargs[arg]})
            else: scatter_kwargs.update({arg: kwargs[arg]})

        ax.streamplot(X_grid[0], X_grid[1], V_grid[0], V_grid[1], **stream_kwargs)

        size = 8 * default_size(adata) if size is None else size
        ax = scatter(adata, layer=layer, color=color, size=size, title=title, ax=ax, zorder=0, **scatter_kwargs)

        savefig_or_show(dpi=dpi, save=save, show=show)
        if not show: return ax
In [26]:
scv.pl.velocity_embedding_stream(adata,
                                 linewidth = 2,
                                 min_mass = 0,
                                 basis = "umap",
                                 color = ['finalclusters','age','condition'],
                                 palette = ['green','red','black','blue','orange','gold','brown','darkgrey'],
                                 alpha=.5,
                                 legend_loc = 'right margin',
                                 legend_fontsize = 20,
                                 title = "RNA Velocity within microglia-tsne",
                                 fontsize = 10,
                                 save = "stream.png",
                                 arrow_color = "black"
                                )
saving figure to file ./figures/scvelo_stream.png
In [22]:
scv.pl.velocity_embedding(adata, basis='umap', legend_loc = 'right margin',color = ['age'], palette = ['lightblue','darkblue'],size = 1, arrow_length=8, arrow_size=2, dpi=300)
In [23]:
scv.pl.velocity_embedding_grid(adata,
                               color=['Meg3'], 
                               layer=['velocity', 'spliced','unspliced'], 
                               arrow_size=5,
                               arrow_length = 10,
                               alpha = 0.2,
                               perc=95,
                               colorbar=True)
In [24]:
scv.tl.rank_velocity_genes(adata, groupby = 'finalclusters')
adata.uns['rank_velocity_genes']['names']
ranking velocity genes
    finished (0:00:00) --> added 
    'rank_velocity_genes', sorted scores by group ids (adata.uns)
Out[24]:
rec.array([('Prc1', 'Kntc1', 'Bex2', 'Abca1', 'Ggt5', 'BC051537', 'Oasl1', 'Clec4a1'),
           ('Nuf2', 'Mybl2', 'Obsl1', 'Pnpla7', 'Ctsf', 'Tmem119', 'Sp100', 'Aoah'),
           ('Sgo2a', 'Brca1', 'Unc5a', 'Gm13710', 'Ecscr', 'Gm15892', 'Gm4951', 'Mndal'),
           ('Kif4', 'Dtl', 'Slitrk3', 'Tcp11l2', 'Olfml2b', 'Trim30a', 'Mx1', 'Ms4a7'),
           ('Hmmr', 'Cenph', 'Syt17', 'Cd300a', 'Gm13710', 'Slc46a3', 'Ifi206', 'Ms4a4a'),
           ('Prr11', 'Gmnn', 'Dsel', 'P2rx4', 'Bmf', 'A3galt2', 'Oasl2', 'Nrap'),
           ('Aurkb', 'Chaf1b', 'Thy1', 'Selenop', 'Ptgs1', 'Selenop', 'Trim30a', 'Fcgr2b'),
           ('Stil', 'Dck', 'Kcna4', 'Ldhb', 'Cmtm8', 'Cd300a', 'Cmpk2', 'F13a1'),
           ('Ankle1', 'Tk1', 'Chd3os', 'Lrrc25', 'Crybb1', 'Lpcat2', 'Samd9l', 'Gm4951'),
           ('Ube2t', 'Dscc1', 'Lhx6', 'Slc46a3', 'Mcm3', 'Il21r', 'Oas1a', 'Pla2g7')],
          dtype=[('1', '<U50'), ('2', '<U50'), ('3', '<U50'), ('4', '<U50'), ('5', '<U50'), ('6', '<U50'), ('8', '<U50'), ('9', '<U50')])
In [25]:
batch1 = bdata[bdata.obs["age"] == "P5",:]
batch1.obs_keys()
Out[25]:
['nCount_RNA',
 'percent.mito',
 'age',
 'condition',
 'finalclusters',
 'sample_description',
 'sex',
 'initial_size_spliced',
 'initial_size_unspliced',
 'initial_size',
 'n_counts']
In [19]:
scv.tl.terminal_states(adata)
scv.pl.scatter(adata, color=[ 'root_cells', 'end_points'])
scvelo.tl.latent_time(data, vkey=velocity, min_likelihood=0.1, min_confidence=0.75,
min_corr_diffusion=None, weight_diffusion=None, root_key=None,
end_key=None, t_max=None, copy=False)
                      
scvelo.tl.velocity_pseudotime(adata, vkey=velocity, groupby=None, groups=None,
root=None, end=None, n_dcs=10, use_velocity_graph=True,
save_diffmap=None, return_model=None, **kwargs)

scv.tl.velocity_confidence(adata)
scv.pl.scatter(adata, color='velocity_confidence', perc=[2,98])
  File "<ipython-input-19-96bdcd23df7e>", line 3
    scvelo.tl.latent_time(data, vkey=’velocity’, min_likelihood=0.1, min_confidence=0.75,
                                              ^
SyntaxError: invalid character in identifier
In [ ]: